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To understand the mechanism of Mott transitions in case of no magnetic influence, superfluid- 
insulator (Mott) transitions in the S = Bose Hubbard model at unit filling are studied on 
the square and triangular lattices, using a variational Monte Carlo method. In trial many-body 
wave functions, we introduce various types of attractive correlation factors between a doubly- 
occupied site (doublon, D) and an empty site (holon, H), which play a central role for Mott 
transitions, in addition to the onsite repulsive (Gutzwiller) factor. By optimizing distance- 
dependent parameters, we study various properties of this type of wave functions. With a hint 
from the Mott transition arising in a completely D-H bound state, we propose an improved 
picture of Mott transitions, by introducing two characteristic length scales, the D-H binding 
length £dh and the minimum D-D exclusion length £dd- Generally, a Mott transition occurs 
when £dh becomes comparable to £dd- In the conductive (superfluid) state, domains of D-H 
pairs overlap with each other (£ d h > £dd); thereby D and H can propagate independently as 
density carriers by successively exchanging the partners. In contrast, intersite repulsive Jastrow 
(D-D and H-H) factors have little importance for the Mott transition. 

KEYWORDS: Mott transition, Bose Hubbard model, two dimensions, superfluid, insulator, variational 
Monte Carlo Method, doublon-holon binding 



1. Introduction 

After early studies of superfluid-insulator or Mott 
transitions in interacting Bose systems, 1 - 1 an experimen- 
tal examples have been realized using ultracold dilute 
gases of bosonic atoms in various optical lattices. 2 ~ 5 - 1 
The essence of these systems is captured by a spinlcss 
Bose Hubbard model with a harmonic confinement po- 
tential. 6 ' 7 ' Aside from this one-body potential, this basic 
model is important on the theoretical ground that a Mott 
transition can be studied without cares of magnetic in- 
fluence, unlike typical fermionic cases. Most researchers 
are certain that, for the Fermi Hubbard models, metal- 
insulator transitions take place at infinitesimal correla- 
tion strength on the hypercubic lattice in any dimen- 
sion. 8 ) Such transitions at U/t = (U: onsite-interaction 
strength, t: hopping integral) require elements other than 
competition between itineracy and bare interaction of 
particles, like magnetic correlation. In contrast, spinless 
Bose Hubbard models bring about Mott transitions at 
moderate finite values of U/t. 9 ~ 12 ^ 

To date, Mott transitions in the S = Bose Hubbard 
model have been studied with various methods. For the 
square lattice of our interest, properties of T c , superfluid 
density, etc. were studied, applying a quantum Monte 
Carlo (QMC) method earlier to small systems 13 ) and 
later to larger systems; 14,15 - 1 a ground-state phase dia- 
gram in the plane of chemical potential and interaction 
strength was constructed using a strong-coupling expan- 
sion. 16 ) These studies estimated the critical interaction 
strength of Mott transitions at U c /t — 16.4-16.7 for the 
particle density of unit filling, n = 1 (n = N/N s with N: 
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particle number, iV s : site number) at T — 0. Thus, the 
existence of a Mott transition has been embodied, but 
the mechanism of the transition is still not clear. 

Variational Monte Carlo (VMC) approaches 17-19 -* are 
very useful to analyze the physics of Mott transitions, be- 
cause one can directly and exactly treat wave functions 
for any values of U/t. As variational approaches to the 
Bose Hubbard model, wave functions with only onsite 
correlation factors, which correspond to the celebrated 
Gutzwiller wave function 20 ) (GWF, ^q) for fermions, 
were studied first. 21,22 ' In contrast to for fermions, 
GWF for bosons is solved analytically without additional 
mean-field-type approximations 23 ) in arbitrary dimen- 
sions, and yield a Brinkman-Rice-type (BR) superfluid- 
insulator transition 24 ' at finite U (= Ubr, see Fig. 2). 
In the insulating side of [/brA, however, all the lattice 
sites are occupied with exactly one particle and the hop- 
ping completely ceases, namely, - > rii=i^]|0) anQl 
the total energy vanishes (E — 0). This result, which 
apparently contradicts experiments and reliable theories 
like the strong-coupling expansion, is caused by an over- 
simplified setup of the wave function, in which the effect 
of density fluctuation should be included. To remedy this 
drawback, it is crucial to add appropriate intersite cor- 
relation factors to GWF. In this line, recent VMC stud- 
ies 25 -* emphasized the importance of an ordinary type of 
long-range Jastrow factor for the transition. 

In this paper, we first show that an attractive correla- 
tion factor between a doubly-occupied site (doublon: D) 
and an empty site (holon: H) plays a leading role to in- 
duce the Mott transition for the present bosonic model. 
In fact, D-H near-neighbor correlations have long been 
studied for Fermi Hubbard models, 26 ~ 29 ) and the present 
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authors have shown that near-neighbor D-H binding fac- 
tors and their analogs are capable of inducing Mott 
transitions for attractive Hubbard models, 30 ) a repulsive 
Hubbard model on an extended square lattice 31 ' 32) and 
on an anisotropic triangular lattice. 33 -* In this mecha- 
nism, in the insulating regime, a doublon and a holon 
as density carriers are confined in the range of near- 
neighbor sites; namely, density fluctuation is localized. 
In this work, we extend the D-H binding correlation fac- 
tor to various long-range types in order to corroborate 
the above conclusion, and propose a renewed picture of 
conduction, which can explain a Mott transition arising 
in a completely D-H bound state. Namely, the Mott tran- 
sition occurs, when a D-H binding length £dh regulated 
by the D-H attractive correlation becomes comparable to 
the minimum D-D (H-H) distance £dd- This corresponds 
to the condition that a D-H pair comes to stop exchang- 
ing the partner with nearby D-H pairs. Furthermore, we 
introduce a repulsive Jastrow factor to check that it is 
of no importance for the Mott transition. These proper- 
ties for bosons are fundamentally common to fermions, 
unless magnetic correlation is explicitly introduced. 34 ) 

This paper is organized as follows: In §2, the model 
and the trial wave functions used in this paper are intro- 
duced. In §3, the properties of the Bose Hubbard model 
as to the Mott transition are studied with a short-range 
D-H binding factor. In §4, we consider the effect of long- 
range D-H attractive correlation factors, and of D-D re- 
pulsive factors. In §5, we discuss a renewed picture of the 
Mott transition by introducing two characteristic length 
scales. In §6, a concise summary is given. In Appendix, 
we briefly note the setup and condition of the VMC cal- 
culations implemented in this paper. 

Parts of the results in this paper have been published 
before. 35 ' 36 ) 

2. Formalism 

After defining the model in §2.1, we introduce trial 
wave functions with various types of D-H attractive cor- 
relation factors and a long-range D-D (and H-H) repul- 
sive factor in §2.2. 

2.1 Bose Hubbard Model 

We consider the S = Bose Hubbard model on the 
square (SQL) and triangular (TAL) lattices with only 
homogeneous nearest-neighbor (NN) hopping: 

H = -t£ + b]h) + ^ ]T n i( n i - !)> (!) 
to) i 

where bj (rij = bjbj) denotes an annihilation (number) 
operator of a boson on the site j. We assume t, U > 0, 
and use t as the energy unit. In this paper, we disregard 
a harmonic confinement potential, not only because we 
allow for the comparison with electron systems in solids, 
but because flat confinement potentials will be expected 
in forthcoming experiments of cold atoms. 37 - 1 We use sys- 
tems of L x L (= N s ) sites with the periodic boundary 
conditions in both x and y directions. The kinetic part, 



H t , in cq. (1) is diagonalized as, 

#t = 5>i4&k, (2) 

k 

with 

_ J — 2t(cosfc x + cosky), SQL 
£k | — 2t [cos k x + cos k y + cos^ + k y )] , TAL 

(3) 

using a Fourier transformation, 

For U/t = 0, since all the particles condense into the 
lowest-energy level k = = (0,0), the ground state is 

$o = Co |0>, (5) 

with the eigenenergy per site being E = —At (— 6t) for 
SQL (TAL). Using eq.(4), the real-space representation 
of <I>o is written as, 

$ cx^6 t 1 4---6j v |0), (6) 
M 

where {r} indicates the sum of all the particle config- 
urations {ri, • • • , tat}, permitting duplicate counting in 
combinations, and we abbreviate 6j. to b\. Note that, 
in contrast to fermionic ground states, the coefficient (a 
permanent of N x N) of every configuration in the sum 
of eq. (6) becomes an identical constant, because all the 
elements in the permanent becomes unity. Accordingly, 
VMC calculations are greatly simplified. 

As mentioned in §1, the ground-state phase diagram 
and some relevant properties of this model including the 
Mott critical value, U c /t, have been studied by reliable 
methods like QMC. Here, we shed light on the properties 
and mechanism as to the Mott transition, taking advan- 
tage of the VMC method. 

2.2 Variational wave functions 

We adopt a many-body variational approach to tackle 
the Mott physics in eq. (1). The simplest trial function is 
a Bose analog of the Gutzwiller wave function 20 ) (GWF), 
*G = g D &o, where g is a variational parameter control- 
ling the number of multiply occupied sites (multiplon: 
M) and D is onsite correlation operator: 

D = \Y, n ^- 1 )- ( 7 ) 

3 

The correct energy expectation value by is given by 
an analog of the Gutzwiller approximation formula, 21, 22 ' 
and a Brinkman- Rice-type transition occurs at ZTbr/ t — 
(I + \/2) 2 z, where z is the number of the NN sites: 
C/brA = 23.31 for SQL, and 34.97 for TAL. How- 
ever, the description of the insulating state by v&q is 
incorrect in that every site is occupied by a single par- 
ticle and density fluctuation is completely suppressed, 
(H t ) — {Hjj) — {Hjj = UD), in the same way as the 
Brinkman-Rice transition in the fermionic cases. 24 ) 
At unit filling, a multiplon and a holon are regarded as 
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positive and negative particle-density "carriers" respec- 
tively in the background of singly-occupied sites of the 
neutral or average particle density. Thus, conductivity 
depends on whether the motion of such density carri- 
ers is free or bound; to describe the Mott transition with 
higher fidelity, it is crucial to add an inter-carrier correla- 
tion factor, especially, a doublon(multiplon)-holon (D-H) 
factor, Tq: 29 1 



DH 



V Q ^ G . 



(8) 



For bosons, multiple site occupation of p (> 0) particles 
is allowed; for U/t = Q, the distribution of the number of 
sites occupied by p particles, P(p), should be Poissonian. 
However, for large values of U/t like in a Mott critical 
region of our interest, P(p) with p > 3 almost completely 
vanishes, as will be discussed later. Then, a multiplon be- 
comes identical to a doublon, and the particle-hole sym- 
metry is restored at n = 1. In this context, we treat 
multiplons and holons symmetrically in Vq, and often 
regard M as D for U £ U c /2. 

In our previous studies for fermions, 29-33 ) we have 
used only near-neighbor D-H binding factors [(1) below], 
which are sufficient to describe Mott transitions. In this 
work, we introduce long-range types of Vq, whose neces- 
sity has long been recognized in an exact-diagonalization 
study in one dimension, 29 ) to study the effect of D-H 
factor more in detail. We itemize Vq used in this paper 
below. 

(1) We extend the above near-neighbor D-H binding 
factors for fermions 38 ) to the Bose Hubbard model. Fol- 
lowing the previous papers, we call this short-range D-H 
binding wave function QWF or tf^' 5 . For SQL, we con- 
sider up to the second (diagonal)-neighbor correlation: 

V Q = V Q (fi,fx') = (l-fxf(l-^f (9) 
where primes indicate diagonal neighbors and, 



Q 



O _ 



E 



+t('), 



(10) 

Here, di and hi are projection operators of multiplon (M) 
and holon (H) on the site i, respectively: 



di\i) 
hi\i) 



1 I*)) 


|*), 



M 

otherwise 



H 

otherwise 



(11) 
(12) 



' runs all the NN (diagonal-neighbor) sites of the site 
i, and |i"isa variational parameter controlling the bind- 
ing of M and H between NN (diagonal-neighbor) sites. 
For TAL, we take account only of the NN M-H correla- 
tion: 



V Q (p) = {\-^ 



(13) 



In short, Q counts the number of M without NN H plus 
that of H without NN M. The range of fi is limited to 
< /i < 1, whereas a subsidiary parameter // sometimes 
becomes negative. 39 - ) The density of M (H) isolated from 



H (M) is reduced by \i; for \i = [p! =) 0, ^p^jj is reduced 
to ^>q, and M can move freely from H. In the other limit, 
H = 1, M (H) cannot appear unless it is accompanied by 
at least one H (M) in the adjacent sites. Such a complete 
D-H bound state is not necessarily insulating, as we will 
see later. In Fig. 1, we show the weight of Vq(/i,0) as a 
function of nearest M(D)-to-H distance r for an interme- 
diate value (fi = 0.7). 




Fig. 1. (Color online) Weight of long-range doublon(D)-holon(H) 
correlation factor in (2) as a function of distance between D and 
the nearest H for several values of the parameter A, eq.(15); (a) 
exponentially decaying type in ^rSj an< ^ P->) power-law decaying 
type in 'J'qjj • For comparison, f(r) of the nearest-neighbor factor 
Vq{0.7, 0) [eq.(9)] in is plotted with dash-dotted lines. 



(2) We introduce long-range D-H binding factors of a 
simple form, which is formally written as, 

Vq{\) = J] { [l - (1 - /(r,)) d 3 ] [1 - (1 - f(r 3 )) h 3 

(14) 

where A is a variational parameter included in f(r), 
which controls the size of an M-H pair. In Vq(\), we 
consider M-H correlation only between each M and its 
nearest H, and vice versa; r (> 1) denotes the distance 
between such M and H in unit of lattice constant, and is 
measured by the stepwise or "Manhattan" metric. This 
fashion suits the spirit of strong-coupling expansion, and 
is not identical with ordinary Jastrow factors, in which all 
the pairs are taken into account. For f(r), on the anal- 
ogy of a fermionic case, 29,34 -* we assume two primitive 
decaying forms, imposing /(l) = 1: 



exp 
1 



1 



(a) exponential 

(b) power 



(15) 



whose behavior is sketched in Fig.l. We write and 
^dh l° r t ne wave functions using the factors of eqs. (15a) 
and (15b), respectively. Compared with the exponential 
form (a) , the power-law form (b) naturally has a tail for 
large r for intermediate to large A c and 1/A p . When A c 
or 1/Ap — > oo, Vq{\) becomes unity, namely, $dh is re- 
duced to ^g- Meanwhile, in the limit of A c or 1/A p — > 0, 
Vq(X) is reduced to Vq((i = 1,/i' = 0), indicating the 
complete D-H binding within the NN sites. Since the 
probability density of M-H pairs with distance r is con- 
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trolled mainly by g for small U/t and by A for intermedi- 
ate and large U/t, the mean M-H distance relates to A c 
or A p for U ^ U c . 

(3) Instead of such an a priori form of f(r) as eq. (15), 
we optimize /(r) in eq. (14) for every r as variational 
parameters, using recently-developed optimization tech- 
niques. 19 ) In this case ['Pq(/)] ) the number of varia- 
tional parameters becomes equivalent to the linear di- 
mension of the system L: g and /(2)-/(L). We represent 
this optimized-/(r) wave function as ^^jh- As we will 
see later, parameters f(r) with large r ( ^ 7), are ac- 
tually unnecessary for any value of U/t, because such 
long-distance D-H pairs rarely appear. Although Vc}(f) 
is naturally better than both Vq{h, //) and Vq(X), we 
can learn much from the comparison with various types 
of TV 

Having considered an attractive part of the long-range 
correlation factor, we should check the effect on a Mott 
transition of a repulsive part, which is known to be im- 
portant for off half filling. 29 > To this end, we consider 
the form, *j = "Pj^dh, with a common repulsive Jas- 
trow factor, 

Vj(a,n) = Y[{[l-(l-r,(r))djd j+r ] 



Thereby, the exponential part of eq. (18) becomes, 



3,r 



x[l- (l-ryM^+r]}, (16) 

in which we take account of the correlation between all 
the M-M and H-H pairs regardless of pair distance, as in 
ordinary Jastrow factors. Here, we adopt a simple power- 
law decaying form of rj{r) as, 



r)(r) = 1 - 



(r > 1), 



(17) 



where a and k are variational parameters with ranges 
< a < 1 and < k < oo. Correspondingly, for Vq in 
$dh, we employ the power-law decaying form in 
eq. (15b). Although the form of eq. (17) may not be the 
best, it is sufficient to ascertain the importance of the 
repulsive Jastrow factor for the Mott transition. 

Incidentally, we compare the above wave functions 
with that studied in related papers by Capello et al.: 25 ^ 



exp 



Vq^o, (18) 



in which they emphasized importance of the long-range 
Jastrow factor. Equation (18) already has an attractive 
D-H factor in the Jastrow (exponential) part, but at each 
the weight of attractive D-H correlation is insepa- 
rably connected with that of repulsive D-D (H-H) corre- 
lation. For U ^ U c /2, since the site occupation number 
is almost restricted to < p < 2, the number operator is 
written as, 



rij — 1 + dj — hj , 



using a doublon operator di in the space of p < 2: 
1 \i), D 



di\i) 







otherwise 



(19) 



(20) 



exp 



o X/ Vi ' j ( didj + _ > - hidj) 



(21) 



which shows the attractive factor is a reciprocal of the 
repulsive factor. Thus, to adjust the D-H binding effect 
independently of the D-D correlation, one is obliged to 
add a redundant Vq term. In this context, have the 
advantage in distinguishing the effect of D-H binding fac- 
tors from that of repulsive ones, although eq. (18) and 
may work similarly. 
We first optimize these trial functions, applying an op- 
timization VMC scheme to systems with up to 1,600 par- 
ticles (L = 40). With the optimal parameters obtained, 
we calculate the expectation values of relevant physical 
quantities, using a conventional VMC method. With this 
procedure, we can obtain accurate variational results in 
most cases, except for statistical errors. In Appendix, we 
briefly explain some details of the VMC calculations car- 
ried out in this paper. 

3. Short-Range Doublon-Holon Factor 

In this section, we study the short-range D-H binding 
wave function QWF, which exhibits typical properties of 
D-H-binding types of wave functions. In §3.1, we study 
the energy of QWF to find out a Mott critical behavior, 
by comparing with extreme cases. In §3.2, we discuss the 
site-occupation number versus U/t, in connection with 
experiments of a quantum gas microscope. In §3.3, the 
existence of Mott transition is corroborated and its prop- 
erties are studied by various quantities. 

3.1 Overall behavior of QWF's energy 

To begin with, we briefly check the optimized energy 
of QWF, eqs. (8)-(13). In Fig. 2(a), the total energies 
per site, E/t, are compared between GWF and QWF of 
Vq(ji, //) for SQL. For small U/t (roughly U < U c /2 ~ 
lOt), the improvement of E/t on GWF is very small, indi- 
cating the D-H correlation plays a minor role for weakly 
interacting conductive states. For U ^ U c /2, however, 
the curve of QWF gradually departs from that of GWF 
and approaches the extreme case of Vq(1,0), in which a 
doublon(s) and a holon(s) rigidly adhere to each other 
in the NN sites. Because such a state strongly suggests 
an insulator, a Mott transition is expected to occur near 
the value where the two curves join (U/t ~ 20). The fact 
that the energy of QWF is considerably lowered from 
that of GWF in this regime indicates the D-H binding 
effect play a prominent part in the Mott physics. The 
transition arising in GWF is a continuous type, because 
E/t vanishes as oc (1 — U /Ubr) 2 in the Bose fluid side. 
In contrast, QWF exhibits a first-order transition, as we 
will see in §3.3. Hence, the mechanisms of the two tran- 
sitions are qualitatively distinct. 

The behavior in the regime of large U/t should obey an 
effective Hamiltonian of eq. (1) in the limit oit/U —¥ 0, 
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which is led by a canonical transformation 40 ) e lS with 

^ = Jj X] {dMbjSj + djtfjbiSi - Sib\bjdj - Sjb^bidi 
(id) 

(22) 

where Si is a projection operator of a singly-occupied 
site: 



1 singly occupied 

\i), \i): otherwise 



(23) 



Then, we have an expression at n — 1 in the space with- 
out D and H as, 



H e tt = 



— ^ SiSj = -2zN a -, 



(24) 



which is drawn as dash-dotted lines in Figs. 2(a) and 
2(b). The energy of QWF well coincides with eq. (24) 
for a wide range of U (> U c ), meaning density fluctua- 
tion is properly introduced in the insulating regime. In 
Fig. 2(b), the same quantity is plotted for TAL. Because 
the behavior is qualitatively identical, henceforth, we ad- 
dress only SQL in most cases. 



U RR /t = 23.31 




n= 1.0 



GWF QWF ft=l L 
• 10 

—a *- 26 

— ♦ — 32 

Square Lattice . 4Q 



10 



20 



30 40 
U m /t = 34.97 




U/t 

Fig. 2. (Color online) Energy expectation values of GWF (open 
symbols) and QWF (solid symbols) as a function of correlation 
strength U/t for (a) the square and (b) the triangular lattices. 
The result of strong-coupling expansion, eq. (24), is drawn by 
the dash-dotted lines. The critical values of Brinkman-Rice-type 
transitions in GWF are indicated by arrows on the upper axes. 
In (a), the data of QWF for Vq(h = 1,/Lt' = 0) are plotted 
with half-solid symbols (expressed as "/j = 1"). The system-size 
dependence is inconspicuous in this scale. 



3.2 Site-occupation number and parity correlation 

In connection with the Mott transition, a change in the 
distribution [P(p)] of site-occupation number p have been 
directly observed in recent experiments of cold bosonic 
atom gases. 41,42 ) Since noninteracting bosons are ran- 
domly distributed to the lattice sites, the ratio of the 
number of sites occupied by p particles (N p ) at U/t = 
should obey a Poisson distribution: 



N, 



N. 



pi 



(25) 



with n — 1 at unit filling. As U/t increases, however, N p 
with large p rapidly decreases to reduce (Hjj) = U(D), 
and for U > U c a number squeezed state is considered 
to be realize, in which most sites are loaded with one 
particle, P(l) ~ 1. The U/t dependence of P(p) has 
been addressed using GWF, 21 ) a mean field approxima- 
tion 43 ) and QMC with analytic confirmations. 44 ) With 
these studies in mind, we discuss the results of $dh- 

Figure 3(a) shows the evolution of P(p) for QWF, as 
U/t is varied. The Poissonian is changed to a shape of a 
symmetric peak centered at p = 1 at a relatively small 
value of U /t, and P(p) with p > 3 almost vanishes; at 
U/t = 12, P(3) is less than P(2)/100 (see Table I). The 
U /t dependence of P(p) depends very slightly on the type 
of D-H correlation factor, as shown in Table I, where the 
values of P(p) are compared among various correlation 




U/t 

Fig. 3. (Color online) Ratio of the number of sites occupied by 
p particles (N P /N B ) calculated with QWF are shown, (a) as a 
function of p for some U/t, and (b) as a function of U/t for small 
p. In (a) the Poisson distribution eq. (25) is plotted with large 
open circles to confirm the VMC outcome for U/t = 0. System- 
size dependence is negligible except for the Mott critical region, 
whose magnification for the case of p = 1 is shown in the inset 
in (b) for several L. 
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Table I. Comparison of probability of site occupation P(p) with 
p = 0-3 among GWF and four types of 'I'dh f° r three values 
of U/t for L = 20. For U/t > 12, P(4) is less than 10 _s in 
every case. In the state column, SF and MI denote 'superfluid' 
and 'Mott insulating', respectively. Incidentally, P(0) and P(2) 
of GWF must completely vanish for L = oo in the MI regime. 



state 



P(0) 



P(l) 



P(2) 



P(3) 



U/t = (Poisson distribution) 



SF 



0.36788 0.36788 0.18394 0.06131 



U/t = 12 



GWF 


SF 


0.14503 


0.71147 


0.14197 


0.00153 


QWF 


SF 


0.12694 


0.74717 


0.12485 


0.00104 


exp. 


SF 


0.13053 


0.73986 


0.12869 


0.00092 


pwr. 


SF 


0.12742 


0.74605 


0.12565 


0.00089 


opt. 


SF 


0.12670 


0.74755 


0.12481 


0.00095 



U/t = 19 



GWF 


SF 


0.04958 


0.90084 


0.04956 


0.00001 


QWF 


SF 


0.02982 


0.94038 


0.02979 


0.00002 


exp. 


MI 


0.01471 


0.97057 


0.01471 


0.00000 


pwr. 


SF 


0.02866 


0.94269 


0.02864 


0.00001 


opt. 


SF 


0.02855 


0.94290 


0.02854 


0.00001 


U/t = 24 


GWF 


MI 


0.00108 


0.99783 


0.00108 


0.00000 


QWF 


MI 


0.01016 


0.97968 


0.01016 


0.00000 


exp. 


MI 


0.01014 


0.97973 


0.01013 


0.00000 


pwr. 


MI 


0.01035 


0.97931 


0.01035 


0.00000 


opt. 


MI 


0.01038 


0.97925 


0.01038 


0.00000 



factors for some values of U/t. Furthermore, the U/t de- 
pendence shown in Fig. 3(b) is quantitatively consistent 
with a result of QMC [Fig. 1(b) in ref. 44]. Thus, in the 
Mott critical regime {U ^ U c /2), we can safely consider 
the problem in the restricted space of p < 2, similarly to 
Fermi systems. 

As seen in Fig. 3(b), the ratio of singly-occupied sites 
(doublon and holon densities) increases (decrease) almost 
linearly with U/t until near the critical point. The vari- 
ation of -P(l) at the Mott critical point is as small as 1% 
in QWF [see the inset of Fig. 3(b)], and local number 
fluctuation remains in some degree even in the Mott in- 
sulating phase. Thus, a number squeezed state is gradu- 
ally approached as U/t increases, and does not distinctly 
characterize a Mott insulating (MI) state. 
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Fig. 4. (Color online) Real-space parity correlation function for 
various values of U/t calculated with long-range D-H binding 
wave function along path shown on abscissa. The arrow in- 

dicates the nearest-neighbor site. The data in the Mott insulating 
regime are indicated by solid symbols. 



Recent development of single atom-single lattice site 
imaging technique (a quantum gas microscope) enables 
one to directly observe the parity (even or odd) of the 
occupied particle number p at each site; 45,46 - 1 the parity 
of the site j is written as, 



1 



-1)"']. 



(26) 



Thus, the correlation function of the parity, 

Cp(j,£) = (pm) - (PMPt) (27) 

becomes a quantity directly measured by experiments. 47 ' 
In Fig. 4, we show C p (j,£) for various values of U/t us- 
ing the long-range ^dh with Vq(f). Because C p (j,£) is 
closely related to the usual density correlation function, 



= (njnt) - (nj)(m), 



(28) 



the magnitude of C p (j,£) rapidly decays, as |r| increases. 
For large U/t, in particular in the insulating regime, 
the parity correlation is almost restricted to the nearest- 
neighbor sites, suggesting insignificance of the long-range 
part of correlation. In this regime, the parity operator 
eq. (26) is written by the projection operators eqs. (12) 
and (20) as, 



(29) 



and rij by eq. (19), so that C p (j, tj is written with N (j, £) 



and the doublon-holon correlation function: 

C w (j,£) = {d j he}-{d j }(he), 



as. 



C p (j,£) = N(j,e) + 4Cjyn{jJ)- 



(30) 



(31) 



Thus, the density and D-H correlation functions are fun- 
damental also in analyzing quantum gas microscope ex- 
periments near the Mott transition. 

3.3 Mott transition in QWF 

Now, we analyze the superfluid-insulator transition in 
QWF more in detail. When we carefully watch the be- 
havior of E/t in a magnified figure (Fig. 5), we notice a 
cusp for each L with L > 20. In each side of the cusp, we 
can find another local minimum (a metastable point) of 
E/t, which is smoothly extrapolated from the curve in 
the other side of the cusp, suggesting a first-order transi- 
tion. This is supported by the existence of a discontinuity 
at the cusp point in the optimized variational parame- 
ters, as shown in Fig. 6. Critical values thus obtained are 
listed in Table III; the system-size dependence of U c /t 
will be discussed in §4.2. A clear sign of a first-order 
transition is not observed for small systems (L < 14 in 
this case), similarly to fermionic systems. 32 ' For TAL, a 
clear discontinuity does not appear even for L = 26, as 
shown in Fig. 7 for /x, and E/t has very broad minimum 
at U ~ U c ; the tendency toward a continuous transition 
is stronger for TAL, which tendency is similar to those 
of frustrated metallic states in fermionic models. 32 ' 

In contrast to GWF, which has g = for U > Ubr, 
QWF has finite g even for U > U c , indicating the exis- 
tence of density fluctuation even in the MI regime. The 
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nearest-neighbor D-H parameter /i exhibits a large dis- 
continuity at U c , and becomes close to 1 for U > U c ; 
the D-H binding becomes firm in the insulating regime. 
Note that, as we will see, the Mott transition is induced 
by the collaboration of the suppression of onsite density 
fluctuation by g and the D-H binding effect by Vq. 

To see the properties of this first-order transition, we 
discuss some relevant quantities. We start with the ex- 
pectation value of onsite correlation operator D, eq. (7). 
Especially at unit filling ((rij) = (nj) 2 =l), this quantity 



-0.3 



-0.35 



Square Lattice 
QWF 

n = 1.0 




U/t 

Fig. 5. (Color online) Magnification of total energy of QWF 
[Fig. 2(a)] near Mott critical points as function of interaction 
strength for several values of L. The arrows indicate the transi- 
tion points (U c /t). Open symbols near U c /t denote the data of 
metastable states. 



0.15 



0.05 




Fig. 6. (Color online) Optimized main variational parameters 
in QWF for square lattice, (a) Gutzwiller parameter and (b) 
nearest-neighbor D-H parameter, near Mott critical values for 
some L. In (a), corresponding values of GWF are added with 
an arrow indicating C/br- Open symbols near U c /t denote the 
values of metastable states. 



per site, d = (D}/N s , coincides with half of the onsite 
density fluctuation or variance, 



N(j,j) 



= (nS ) - (nj) = 2d (n = 1). 



(32) 
(33) 



Substituting eq. 



"3 1 V3I 

(19) in eq. (33), we have a relation, 
d = d, valid for U ii U c /2, where d is the doublon den- 
sity. Thus, d is a key indicator of Mott transitions like 
for fermionic systems, so that we regard d as an order 
parameter of Mott transitions also for bosons. In Fig. 8, 
we show the U /t dependence of d for GWF and QWF, 
and its magnification near U c in the inset. For each value 
of L, d of QWF decreases with U/t and exhibits a dis- 
continuity at each U c /t (L > 20). Conversely, the kinetic 
energy E t = (H t ) increases with a discontinuity at U c 
(not shown). This energetics meets the ordinary criterion 
of Mott transitions: The energy is stabilized by lowering 
Ejj at the cost of E t for U > U c . The doublon density 
d for QWF remains finite in the insulating regime, be- 
cause local density fluctuation is permitted and transient 
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7. (Color online) Optimized nearest-neighbor D-H parame- 
ter near Mott critical values in QWF for triangular lattice. Four 
system sizes (L) are compared. 




U/t 

Fig. 8. (Color online) Comparison of expectation values of on- 
site correlation operator D [eq. (7)] per siteon the square lattice 
among GWF, QWF and completely D-H bound state (dis- 
cussed later). This quantity is substantially doublon density d for 
U/t <; 10, and is half the onsite density fluctuation a\. The ar- 
row on the curve of 1 indicates the Mott critical point. The 
inset shows the magnification near the Mott critical points for 
GWF and QWF. The symbols in the inset are common to those 
in the main panel. 
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D-H pairs exist within neighboring sites [see also P(2) in 
Table I]. This behavior contrasts with that of GWF, for 
which d vanishes for U > Ubr- Note that this quantity 
will be measured as a function of U / t by experiments like 
quantum gas microscopes. 




0.05- 



Fig. 9. (Color online) (a) Condensate fraction as a function of 
U /t for some values of L, calculated with three functions, GWF, 



QWF, and completely D-H bound function, * 



Q 



(discussed 



later). The arrow on >I/q indicates the Mott critical point, (b) 
Magnification of the Mott critical region for GWF and QWF. 



Next, we consider U/t dependence of the condensate 
fraction, 



po = n(0)/N s 



(34) 



with n(0) being the k 
distribution function: 

n(k) = (blb k ) 



clement of the momentum 



1 



(35) 



In Fig. 9, we plot po calculated with GWF and QWF. As 
U/t increases, the condensate fraction diminishes from 
the value of free bosons po = 1, and vanishes at U = U c 
as the order of n/N s in the non-supcrfluid states for 
all the wave functions. Thus, superfluidity vanishes at 
U = U c . As shown in Fig. 9(b), the discontinuities ap- 
pear for QWF with L > 20 in accordance with other 
quantities. Recently, the condensate fraction has been 
actually observed by experiments of cold atoms. 48 ' 49 ' 

Here, we discuss k ^ elements of n(k), because it is 
a directly-observed quantity by cold-atom experiments. 
Figure 10 shows n(k) obtained by the VMC calculations 
under some conditions. Because GWF gives the result 
identical with that of the Gutzwiller approximation in 



fermionic cases, 23 * 1 n(k) is constant for k ^ in both su- 
perfluid and insulating phases. On the other hand, wave 
functions with D-H factors yields dispersive n(k), ow- 
ing to the effect of density fluctuation. In Fig. 10, we 
draw n(k) of the best D-H binding wave function in this 
study It is marked that the difference is small be- 

tween the superfluid (U/t = 20) and insulating states 
(U/t = 24). As k leaves the r(= 0) point, n(k) de- 
creases in any direction, but its decrement is larger in 
the r — > (it, 7r) direction than in the T —> (tt, 0) direction. 
This direction-dependent dispersion near the Mott criti- 
cal point is actually observed as a cross-like intensity in 
absorption images of cold atoms, for example in Figs. 2f 
and 2g of ref. 2. This topic was previously argued by the 
perturbative correction to the Gutzwiller solution. 50 ) 

Finally, to corroborate a superfluid-insulator transi- 
tion at unit filling, we compare particle-density (n) de- 
pendence of various quantities in the vicinity of n = l 51 - 1 
between U <U C and U > U c . Let us start with the con- 
densate fraction po, which is plotted in Fig. 11(a). For 
U = 16t (< U c ), po becomes minimum at n = 1, but pre- 
serves a finite magnitude, whereas for U/t = 24 (> U c /t), 
po decreases as n approaches 1, and almost vanishes at 
n — 1. The sign of second derivative d 2 po/dn 2 becomes 
different between the two cases near unit filling. Thus, 
a MI state appears only at n = 1 and U > U c . This 
is supported by the similar behavior of the onsite den- 
sity fluctuation ct 2 , eq. (32), as shown in the inset of 
Fig. 11(b). The order parameter of Mott transitions, d, 
depicted in the main panel of Fig. 11(b) monotonically 
increases for U < U c , whereas it sharply decreases at 
n = 1 for U > U c , suggesting a Mott transition. Figure 
11(c) represents the chemical potential, 

c(N-^j=E(N)-E(N-l) = ^- (36) 

as a function of n. Although for U < U c , Q/t is always 
a smooth function of n, for U > U c , it has a discontinu- 
ity corresponding to the Mott gap, broadly estimated as 
A/t = 8.2 for U/t = 24 on SQL, and 17.6 for U/t = 40 
on TAL (figure not shown). This is a direct evidence of 
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Fig. 10. (Color online) The momentum distribution function, 
eq. (35), are compared between GWF (^g) an d optimizing wave 
function (^q 1 ^), and between superfluid (SF) and Mott insulat- 
ing (MI) phases. The value of the coherent wave number k=0 is 
omitted. Two system sizes are used. The behavior of n(k) in the 
other D-H binding wave functions is similar to that of ^n^j. 
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Fig. 11. (Color online) Particle-density dependence of various quantities are compared between at U/t = 16 [a superfluid (SF) state 
for n = 1] , and at U/t = 24 [a Mott insulating (MI) state for n = 1] , calculated with QWF on the square lattices of some system sizes 
L X L. (a) Condensate fraction given by n(0)/N a . (b) The main panel shows the substantial doublon density, (D)/N a , and the inset 
the onsite density fluctuation ct^ = N(j,j) defined by eq. (32). For n ^ 1, the relation = 2d [eqs. (32) and (33)] does not hold, (c) 
Chemical potential, eq. (36), actually estimated from the finite difference of E/t with respect to n. The arrow indicates the Mott gap 
for U/t = 24 at n = 1. (d) Optimized nearest-neighbor D-H binding parameter. 

the Mott transition. 

To recognize the importance of the D-H binding effect 
on the Mott transition, we compare, in Fig. ff (d), the 
optimized D-H binding parameter fi in QWF between 
the two phases. In the two phases, \x is symmetric with 
respect to n = 1, and has a maximum at n = 1, but 
the magnitude is distinct. For U = 16t, \i slowly varies 
with n and is still small at n = I, whereas for U = 24t, 
jj, anomalously increases as n approaches 1, and almost 
reaches 1 as L increases. Thus, the D-H binding effect is 
significantly enhanced in the very vicinity of the insulat- 
ing state. 

4. Long-Range Correlation Factors 

In this section, we discuss the effect of long-range cor- 
relation factors, which we have disregarded in the preced- 
ing section. In §4.1 and §4.2, we focus on the properties 
of D-H attractive factors. In §4.3 we study the effect of 
additional D-D and H-H repulsive factors. 



4--1 Optimized D-H attractive factors 

To begin with, we compare the minimized energy 
among QWF and the wave functions with three long- 
range D-H attractive factors, exponentially decaying 
(fpg), power-law decaying (^^jh)' an d completely op- 
timizing (^qjj) types, introduced in §2.2. As will be dis- 
cussed in §4.2, a Mott transition occurs in each wave 
function near U c of QWF (~ 20t). For sufficiently large 



and small U, compared with U c , the total energies of 
the four functions are close to one another. As shown 
in Fig. 12, some difference appears in the Mott criti- 
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Fig. 12. (Color online) The total energies are compared near the 
Mott critical points among the wave functions with four kinds of 
D-H attractive factors: a short-ranged type (QWF) and three 
long-range types (exponentially- and power-law-decaying and 
completely-optimizing types). For clarity, only the data of two 
system sizes are plotted. The Mott critical point for each func- 
tion (L = 20, and 14 for 'i'^jj) are indicated by arrows. The 
value of strong-coupling expansion (—8t/U) is added. 
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cal region. Of course, E opt (E for *° P h) is always the 
lowest. For U/t <^ 21, the energies of three long-range 
D-H wave functions are broadly analogous and some- 
what improved over _gQ WF for each system size. On the 
other hand, for U/t ^ 20, E e * p becomes clearly higher 
than E/t of the other D-H wave functions; E pwr is also 
slightly higher than E I ^ WF . Now, we are aware that the 
frequently- used QWF yields a relatively good result, es- 
pecially for U < U c , despite its simplicity. We analyze 
these features of E/t by comparing the optimized forms 
of correlation factors f(r) among the D-H wave func- 
tions, in the following. 




r U/t 



Fig. 13. (Color online) The appearance probability of a doublon 
the distance from which to the nearest holon is r is plotted, (a) 
for the completely optimizing wave function 'SPjjH {Uc/t ~ 21.3), 
versus r and for some values of U/t, and (b) for the four D-H 
wave functions versus U/t for r = 1 and 2. The threshold value 
Pl) ln = 0.0004 is indicated by a dash-dotted line in both panels. 

First, we discuss the effective range of f(r). We repre- 
sent the distance from a doublon (multiplon) to its near- 
est holon simply by r here, and the probability that a site 
is occupied by a doublon of r by pT>(r), which satisfies 
the relation, 

L 

£>(r) = d. (37) 
r=l 

In Fig. 13(a), we plot pr>(r) of vfr^ in a logarithmic scale. 
Because pjj(r) is a monotonically decreasing function of 
r, it is convenient to define a threshold value p™ m , which 
broadly gives the effective range of f(r) by, 

Mr) > Prf • (38) 

The contribution from pv(r) < pf^ m should be negli- 
gible, namely, the corresponding particle configurations 
should appear very rarely. Thus, the optimized values 
of f(r) for such distant r are unreliable and insignifi- 
cant. Here, we set pg in = 100/250,000 on the basis of 
accuracy in the VMC calculations. In a weakly inter- 
acting regime (U/t Ss 12), the weight of pr>(r) concen- 
trates in r ^ 4, and virtually vanishes (pv(r) < Pd'") f° r 
larger r, because multiplons and holons are crowded. In 
the superfluid regime, pv(r) comes to decrease slowly as 
U/t increases, and extend the valid range to r 6 near 
U c /t. Meanwhile, in the insulating regime, pv(r) imme- 



Table II. Comparison of appearance probability and weight of D- 
H factors for r = 1-3 among four wave functions in insulating 
regime (U/t = 22) for L = 20. 





Pd(1) 


Pd(2) 


/(2) 


Pd(3) 


/(3) 


E/t 




xl0~ 2 


xlO- 4 




xl0~ 5 






optim. 


1.18 


2.54 


0.263 


2.57 


0.090 


-0.3282 


power 


1.17 


2.08 


0.247 


3.80 


0.109 


-0.3279 


exp. 


1.15 


1.75 


0.236 


0.84 


0.056 


-0.3272 


QWF 


1.05 


0.052 


0.050 


0.64 


0.050 


-0.3246 



diately decays to have substantial weight only on r = 1 
or at most 2. This contrastive feature of the effective 
range clearly implies that the nature of D-H binding ef- 
fect changes at U c . This effective range of f(r) is closely 
related to the D-H binding length, £dh, introduced in §5. 
Making a similar analysis with the condition (38) for the 
other wave functions, we determine the effective range of 
f(r) for various values of U/t for each function. 




Fig. 14. (Color online) The optimized D-H attractive factors are 
compared among four wave functions versus nearest D-to-H dis- 
tance r. Except for ^Q^j, the data for U/t = 12, 16, 19 and 
21 are shown for clarity. The range of r is restricted according 
to the condition (38). The data of U/t = 19 for *p X P and of 
U/t = 21 except for ^q 1 ^ are in the insulating phase. For QWF, 
we disregard the contribution of //. 

Figure 14 compares the behavior of optimized f(r) 
among the four wave functions within the effective range 
thus determined. In the conductive regime, f(r) in ^q 1 ^ 
rapidly decreases for r & 3, but becomes almost constant 
for r 3. Namely, the D-H binding is effective only for 
r & 3, and a doublon is released from the bondage of 
holons for r ^ 3. Because f(r) of QWF is constant for 
r > 3, the behavior for large r is analogous to f(r) in 
^p 1 ^. In contrast, f(r) in and especially in 

continues decreasing to zero as r increases, namely, the 
effective range of D-H binding is too long, compared with 
f(r) in vPqh- This is a cause of the unexpected good (un- 
satisfactory) result of QWF (f^g). In the MI regime, 
f(r) in ^y^i rapidly decays with r, and the effective range 
is limited to at most r — 2, as mentioned [Fig. 13(a)]. 
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Regarding energy improvement, the appearance proba- 
bility of nearest-neighbor D-H pairs seems primarily im- 
portant; as listed in Table II, pd(1) is more than 10% 
smaller in QWF than in the long-range wave functions. In 
addition, the role of /(2) is not negligible. In Fig. 13(b), 
Pr>{ r ) for r — 1 and 2 are compared in a logarithmic 
scale. For U > U c , the three long-range wave functions 
have similar values of /Qd(2), whereas the values of QWF 
are roughly two orders of magnitude smaller. This is di- 
rectly reflected in /(2), as shown in Table II. No such 
great differences in pv(r) and f(r) can be seen for r > 3 
among the four wave functions [see also Fig. 15(b)]. 

4-2 Mott transitions in D-H attractive factors 




10 u/ 20 30 




I 1 , 1 , 1 , 1 , 1 

16 18 20 22 24 

U/t 



Fig. 15. (Color online) (a) The optimized parameters controlling 
D-H correlation length, 1/A C in and A p in are plotted 

for L = 10-20. (b) The weight of D-H attractive factors f(r) for 
r = 3 is compared near the Mott critical points among the four 
D-H binding wave functions. In both panels, the critical points 
are indicated by vertical arrows, in case they are manifest. The 
thick arrows for ^q 1 ^ in (b) means large ambiguity in U c /t. 

In this subsection, we consider other properties of long- 
range wave functions as to the Mott transition. 

We start with the parameters which controls the range 
of D-H attractive correlation, namely A c in anc ^ -V 
in (Fig. 1). Their optimized values are shown in 

Fig. 15(a). Both A and 1/A P abruptly increases at U/t = 
17-20, namely, the D-H correlation range becomes short. 
In particular, clear jumps exist for L > 14 in 1/A C and 
L = 20 in A p , as indicated by arrows. In Fig. 15(b), 
the optimized D-H correlation weight f(r) for r = 3 is 
magnified near the critical points. Both 'f^j' and ^j^ 1 



exhibit critical behavior at the same U c /t as for A; Jumps 
in /(3) exist also for Vt^m with L > 20. 




u/t u/t 



Fig. 16. (Color online) (a) Variance of onsite density correlation 
function (or 2d) and (b) condensate fraction, both calculated 
using wave functions with three kinds of long-range D-H factors, 
*dh an( ^ ^dh f° r a f ew va lues of L. The arrows shows 
the Mott critical points. 

Next, in Fig. 16, we show two quantities characteriz- 
ing Mott transitions, er [eq. (32)] or 2d [eq. (33)], and 
po [eq. (34)]. Again, both exhibit anomalous behavior at 
U c /t determined above by the anomalies of variational 
parameters. The behavior is similar to that of QWF 
(Figs. 8 and 9). Thus, a first-order Mott transition can be 
described by a wide class of D-H binding wave functions 
if sufficiently large systems are considered. 

Details of the critical behavior as well as the value 
of Mott critical point are different for different wave 
functions in some degree. As clear cusps appear in E c * p 
(Fig. 12), ^J^ji exhibits sharp first-order critical behavior 
even for a small system of L = 14. Near U c , two energy 
minima are distinguished, and a hysteresis is confirmed in 
the E/t-U/t plane. A hysteresis is also observed for 
widely near U c /t, and the energy difference between the 
two phases is very small. Thus, the control of optimiza- 
tion process becomes considerably difficult in the critical 
regime, so that the estimation of U c /t for is less 

accurate than for the others. This is probably caused by 
the redundancy of parameters. In comparison, be- 
haves mildly near U c and we have not detect a manifest 
hysteresis even for L = 20. 

Finally, we discuss the value of the Mott critical point. 
In Table III, the critical values determined by the D- 
H binding wave functions in this work are listed with 
those of other studies. Among the three long-range V&dh 
treated here, U c /t is mutually different to some extent. 
On the basis of E/t in Fig. 12, this difference is con- 
sidered to stem from the propriety of $dh in describing 
the superfluid state near U c /t. Because E/t of the three 
$dh's is similar to one another for U > U c , a lower 
E/t or a better ^dh in the conductive side of U c yields 
a larger U c /t. Furthermore, the critical value U c /t, es- 
timated by any $dh including QWF (§3.3), steadily in- 
creases as L increases. This is mainly because the system- 
size dependence of E/t is considerably large in the insu- 
lating side of U c /t, but small in the conductive side, as in 
Figs. 5 and 12. These aspects of the critical value are con- 
trary to reliable estimates by QMC and a strong-coupling 
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Table III. The Mott critical values U c /t are compared among the 
wave functions treated in this paper for some system sizes. Sec- 
tions with hyphens indicate the cases in which clear behavior of 
first-order transition is not observed. For empty sections, we have 
not performed calculations in this study. The last low shows the 
result, when the repulsive Jastrow factor is included (see §4.3). In 
the lower lines, the values obtained in other studies are entered. 



\I/ (or method) 


Uc/t (n = 1, square lattice) 


L -> 


10 


14 


20 


26 


32 


40 


QWF (*q) 
optim. (ffgg) 
exp. 

power (^dh) 




17.2 


20.5 
~ 21.3 
17.65 
19.85 


21.5 
~ 22.7 


22.5 


24.1 


repulsive (^j) 






19.85 








GWF (C/brA) 


23.31- ■ • (refs. 21 and 22) 


VMC (Jastrow) 


~ 20.6 (ref. 25) 


Strong coupling 


16.7 (ref. 16) 


Recent QMC 


16.25 (ref. 14), 16.7 (ref. 15) 



phases, and comparable to the statistical errors, which 
are large for 'I'j. 53 ) 




Fig. 17. (Color online) (a) The variance of onsite density correla- 
tion function (or 2d) and (b) condensate fraction, are compared 
between Wgg' (D-H) and *P wr (+D-D) for three values of L. 
The arrows indicate U c /t for L = 20. 



Table IV. Comparison of total energy per site E/t near U c /t 
among GWF and wave functions with power-law decaying type 
D-H and D-D correlation factors. The digits in the brackets in- 
dicate the ratios with respect to the values of GWF+D-H. 



u/t 


GWF 


GWF+D-H 


GWF+D-H+D-D 


L = 14 


19 


-0.11712 (0.301) 


-0.38945 


-0.39156 (1.005) 


20 


-0.07435 (0.200) 


-0.37103 


-0.37561 (1.012) 


21 


-0.04314 (0.121) 


-0.35523 


-0.35680 (1.004) 


L = 20 (U c /t = 19.85) 


19 


-0.11042 (0.293) 


-0.37742 


-0.37858 (1.003) 


20 


-0.06750 (0.191) 


-0.35410 


-0.35405 (1.000) 


21 


-0.03612 (0.106) 


-0.34034 


-0.34138 (1.003) 



expansion, U c /t = 16-17. To obtain a more accurate crit- 
ical value, some factor overlooked in the present class of 
wave functions may be taken into account. As a possibil- 
ity, we take up D-D repulsive correlations in §4.3. 

4-3 Effect of repulsive Jastrow factors 

So far, we have disregarded the effect of intersite re- 
pulsive factors, because it is known long-range repulsive 
factors reduce the energy only slightly for fermions at 
half filling. 29 ' 34 ^ Here, we check that the Bose Hubbard 
model also has this property, and a repulsive factor is in- 
sufficient to improve the system-size dependence of U c /t. 
For simplicity, we employ the power-law-decaying type 
both for attractive [eqs. (14) and (15b)] and repulsive 
[eqs. (16) and (17)] correlation factors mentioned in §2. 

First, we argue improvement in energy by D-H attrac- 
tive and D-D repulsive factors over GWF. Table IV lists 
the numerical values of E/t for the three wave functions 
and two system sizes. As U approaches Ubr (= 23.314), 
E GWF steadily increases toward zero. By introducing the 
D-H binding factor, E is significantly improved on E GWF 
(70-90%), and comes to increase slowly as U/t increases, 
like the case of QWF in Fig. 2(a). On the other hand, 
when we add the D-D and H-H repulsive correlation to 
f dh, the improvement of E/t on f dh is very slight 
(mostly less than 1%) in both conductive and insulating 



Although the energy reduction is slight, we should con- 
firm whether the D-D repulsive factor Vj improves phys- 
ical quantities relevant to the Mott transition. Since Vj 
tends to lengthen the inter-doublon (D-D) distance in the 
insulating side (not shown) , we have expected that ^ wr 
lowers U c /t. In Figs. 17(a) and Figs. 17(b), we compare 
(T and pa, respectively, between and • Against 
our anticipation, the differences in both quantities are 
negligibly small, and a meaningful shift of U c /t cannot 
be observed (Table III). 

Here, we have focused on the power-law decaying case. 
For fermions, it is found that various repulsive factors 
somewhat lower the values of U c /t, but the situation in 
system-size dependence does not change. 34 - 1 Hence, we 
conclude that the effect of repulsive Jastrow factor is not 
significant, as far as the Mott transition is concerned. We 
will take no account of this factor in the following. 

5. Renewed Picture of Mott Transition 

In a previous paper for electrons, 32 ) we affirmed a sim- 
ple mechanism of Mott transitions, in which the bind- 
ing of a doublon (plus charge carrier) to a holon (minus 
charge carrier) in the insulating regime and the release 
from binding in the conductive regime are the essence of 
the Mott transition. In this section, we extend this pic- 
ture to comprehend a case of completely D-H binding. 

First, we briefly mention the properties of the com- 
pletely D-H bound state in NN sites, ^q -1 , namely QWF 
with /i(= 1) and //(= 0) in eq. (9). Because, in ^q -1 , 
a doublon must be accompanied by at least one holon 
in its four NN sites, we tend to regard it as insulating 
for any value of U/t. In fact, however, ^q -1 is a super- 
fluid state for small U /t, and exhibit a Mott transition at 
U/t = 4.55. A sign of a Mott transition can be recognized 
in the cusp behavior of E/t at U/t ~ 4.55 in Fig. 2(a), 
and especially in sudden drop of the condensate fraction 
in Fig. 9(a). In Fig. 18, we show the optimized value of 
g and multiplon density d, besides the above two quanti- 
ties. The discontinuous decreases of the order parameter 
d, and the sudden vanishing of po corroborate the Mott 
transition at U c /t — 4.55. 
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This Mott transition in ^g -1 , especially in the conduc- 
tive side, cannot be understood by a simple release from 
the D-H binding and mutually independent movement of 
the two kinds of carriers, because a doublon is in contact 
with a holon even in the conductive regime. Instead of 
this simple point of view, we propose an improved picture 
extensively applicable to Mott transitions. To this end, 
it is convenient to introduce two relevant length scales, 
D-H binding length £dn, and minimum D-D (H-H) dis- 
tance £dd, which are loosely defined as follows: The dis- 
tances of most of the nearest D-H pairs are smaller than 
£dh; hi other words, a doublon is seldom distant from a 
holon beyond £dh- £dd is a D-D (H-H) exclusion distance, 
namely, the inter-doublon (inter-holon) distances within 
which two doublons (or holons) are mutually almost in- 
accessible. In general, £dh as well as £dd depends on U/t. 

Then, we postulate that an attractive correlation fac- 
tor Vq yields D-H pairs of a binding length £dh accord- 
ing to U/t. An improved picture of Mott transition is 
schematically represented in Fig. 19. In the insulating 
phase, the relation £dh < £dd holds, indicating that the 
domains of D-H pairs do not usually overlap, at least, not 
in sequence. As a result, most D-H pairs are isolated and 
a doublon and a holon are confined within £dh, resulting 
in only local density fluctuation. To this point, the pic- 
ture is basically identical with one previously proposed. 
In the conductive phase (U < U c ), £dh becomes longer 
than £dd, indicating the domains of D-H pairs overlap 
with one another. Then, a doublon in a D-H pair can 
exchange a partner holon with a holon in an adjacent D- 
H pair. Consequently, a doublon and a holon can move 
independently as carriers by exchanging the partner, as 
shown in long arrows in Fig. 19. It follows that, as the 
value of U/t is varied, a Mott transition takes place when 
£dh becomes equivalent to £dd, which (correctly i^d) is 
roughly l/y/d, and is expected to be a monotonically 
increasing function of U/t. 

To justify the above picture, we have to appropriately 
estimate £dh and £dd, and confirm that the Mott tran- 




U/t 



Fig. 18. (Color online) Four quantities calculated with the com- 
pletely D-H bound state 1 are plotted as a function of U/t; 
namely, total energy per site Eft (right axis), half of the op- 
timized Gutzwiller parameter g, multiplon density d and con- 
densate fraction po- The Mott critical point U c /t is indicated 
by a shadowed dash-dotted line. The system-size dependence is 
negligible. 
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Fig. 19. (Color online) Schematic figure of improved picture of 
Mott transitions. Small solid and open circles denote doublons 
and holons, respectively. The D-D distance is closely related to 
§dd> shown by "dd". Each pale circle indicates the domain of 
a D-H pair, whose diameter roughly corresponds to ^j^i shown 
by "dh". Typical configurations are drawn for the two phases. 
In each pair domain, the existence of at least one doublon and 
one holon is required, but excess doublons and holons can move 
around independently slipping out of the original domain, as 
indicated by long arrows for the conductive phase. 



sitions take place when these two length scales intersect 
each other at U c /t. After checking various cases, we have 
found that the following formulae work well and are phys- 
ically natural to the above point of view, 



£dh 
£dd 



Kdh 



Kdd 



Cdh, 



Odd, 



(39) 
(40) 



where £dh (£dd) denotes the average of the nearest D-to- 
H and H-to-D (D-to-D and H-to-H) distance, and a a the 
standard deviation of £a, with A being an index of "dh" 
and "dd": 



\ 



1 



M ^ 

m— 1 



(41) 



Here, m runs over all the doublons and holons in all 
the measured samples, M indicates the total number of 
doublons and holons in all the measured samples, and 
indicates the nearest A distance for the m-th doublon (or 
holon). The addition and subtraction of a in eqs. (39) 
and (40) represent a U/t- dependent "softness" of the D- 
H binding and D-D repulsive correlations, respectively. 

First, we discuss a special case of the tightly D-H 
bound state ^q -1 , shown in Fig. 20(a). Because ^dh is 
restricted to 1 in *q _1 , Odh = and £dh = 1 hold, irre- 
spective of the value of U/t. On the other hand, because 
the doublon density decreases as U/t increases, both Idd 
and Odd monotonically increase qualitatively in a man- 
ner similar to those of QWF shown in Fig. 21, and conse- 
quently, £dd also monotonically increases with a discon- 
tinuity at U c /t. As shown in Fig. 20(a), £dd intersects £dh 
in this discontinuity; this special case is consistent with 
the above picture. 

Next, we consider more ordinary D-H binding wave 
functions: QWF, and *dh- In Fi S- 21 > we 
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U/t U/t XJm.lt 

Fig. 20. (Color online) The D-H binding length fdh and mini- 
mum D-D distance £ dd defined by eqs. (39) and (40) are plotted 
versus U/t for (a) the completely D-H bound state, and (b) the 
Gutzwiller wave function, for some values of L. 




t 



Fig. 21. (Color online) (a) The average nearest D-H and D-D 
distances and (b) their standard deviations are plotted for the 
short-range D-H binding wave function QWF as a function of 
U/t. The Mott critical values for L = 20 and 26 are indicated 
by arrows in (a). The scales of coordinates in (a) and (b) are 
common. 

show the average nearest D-H and D-D distances and 
their standard deviations for QWF. The behavior of £'s 
and of c's for the three long-range D-H binding wave 
functions is basically the same. As U /t increases, £dd 
rapidly increases, exhibits a discontinuity at U c /t, and 
further increases for U > U c with a considerable system- 
size dependence. Meanwhile, i^h at first gradually in- 
creases, and has a maximum in the vicinity of U c /t, 
then decreases discontinuously at U c /t, and converges 
to ^dh = 1 for U > U c . The magnitude of Odh becomes 
large near U c /t in the conductive side, but tends to van- 
ish for U > U c . Substituting such I and a for eqs. (39) 
and (40) , we obtain £dh and £dd for each of the four D- 
H binding states, and plot them in Fig. 22. As expected 
from Fig. 21, £dd is a monotonically increasing function 




10 20 30 10 20 

U/t U/t 



Fig. 22. (Color online) The D-H binding length and minimum 
D-D distance £dd defined by eqs. (39) and (40) are plotted versus 
U/t for (a) the short-range D-H binding state QWF, (b) the 
exponentially-decaying D-H binding state ^lth, (c) the power- 
law-decaying D-H binding state , (d) the optimized D-H 
binding state ^p 1 ^, for some values of L. An arrow in each panel 
denotes the Mott critical point for the largest system size in each 
state. 

of U jt with a discontinuity at U c , whereas £dh has a max- 
imum just below U c and becomes a decreasing function of 
U/t with a discontinuous drop at U c . Consequently, the 
magnitude of the two length scales is reversed suddenly 
at U c , and the relation £dd > £dh holds for U > U c for 
each L of every wave function. Thus, we conclude that 
a picture explained in Fig. 19 seems appropriate to the 
D-H binding mechanism of Mott transitions. 

Incidentally, the above point of view is out of focus for 
a non-D-H binding mechanism like Brinkman-Rice-type 
transitions. As depicted in Fig. 20(b), for GWF, £dh is 
always longer than £dd, and does not show a tendency to 
decrease, up to the Brinkman-Rice point. The value of 
£ at ?7brA is proportional to L, which means that ^>q 
entirely lacks intersite correlations. 

6. Summary 

In this paper, we have studied the spinless Bose Hub- 
bard models at unit filling on the square and triangular 
lattices, using a variational Monte Carlo scheme. Our 
primary aim is to grasp fundamentals of the Mott tran- 
sition without influence of the spin degree of freedom. 
In the trial wave functions, we allow for various types of 
doublon-holon attractive factors and a doublon-doublon 
(holon-holon) repulsive factor in addition to the onsitc 
repulsive (Gutzwiller) factor. We itemize the main re- 
sults below. 

(1) Because multiply occupied sites with more than 
two particles almost vanish for U ^ U c /2 for all the wave 
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functions we check, we can consider the Mott physics 
near U c /t (~ 20 for the square lattice) only with a doubly 
occupied site (doublon, D), an empty site (holon, H) and 
a singly occupied site, like fermionic cases. 

(2) Wave functions with appropriate D-H attractive 
correlation factors exhibit first-order Mott transitions. 
Unlike the Brinkman-Rice-type transition arising in the 
Gutzwiller wave function, these transitions have density 
fluctuation in the Mott insulator side. By density fluctu- 
ation, anisotropy is introduced into the momentum dis- 
tribution function, which is consistent with absorption 
images observed in cold atom experiments. 2 ' The total 
energy for U > U c coincides well with the results of the 
strong coupling expansion (E/t cx —t/U). 

(3) In the conductive (superfluid) state, the optimized 
value of D-H attractive correlation weight f(r) [eq. (14)], 
with r being the interparticle distance, rapidly decreases 
with r for r ^ 3, and is almost constant for larger r. 
Thereby, we found it reasonable that the wave function 
with short-range D-H factors (QWF), in which f(r) for 
r > 3 is constant, gives unexpectedly good results. On 
the other hand, in the Mott insulating state, f(r) almost 
vanishes for r > 3, suggesting that a doublon and a holon 
confine each other in near-neighbor (r < 2) sites as a 
D-H pair. This will be confirmed by recently developed 
quantum gas microscope experiments. 45 ' 46 - 1 

(4) The critical behavior in the triangular lattice is 
more continuous-transition-likc than in the square lat- 
tice. A similar phenomenon emerges in paramagnetic 
states of electronic systems, 32 ' and was ascribed to the 
frustration of spins. However, it is probable that its origin 
is closely related to the connectivity of lattices. 

(5) We have improved the D-H binding picture of Mott 
transitions, by introducing two characteristic length 
scales, D-H binding length £dh, which broadly represents 
the size of a D-H pair, and minimum D-D distance £dd- 
We appropriately determined £dh and £dd- In the conduc- 
tive state (£dh > £dd), because the domains of D-H pairs 
mutually overlap, density carriers (D and H) can move 
independently of the partners of the D-H pairs released 
from the binding. In the insulating state (£dh < £dd) ; be- 
cause most domains of D-H pairs are detached from one 
another, density fluctuation is localized within the do- 
main of £dh- The Mott transition takes place, when the 
relation £dh = £dd is satisfied. 

(6) By adding a D-D (and H-H) repulsive factor, the 
variational energy is improved slightly, especially in the 
insulating regime. However, we have not recognized qual- 
itative influences thereof on the Mott transition. 

Since the properties of the Mott transition studied in 
this paper basically coincide with those for the electron 
systems, 32 ' the renewed picture of Mott transitions given 
in §5 can be applied to electron systems. 34 ' An important 
remaining problem is to determine the Mott critical value 
U c /t more accurately. This problem is closely connected 
with the great system-size dependence of the present trial 
wave functions around U c /t. To remedy it, we may need 
more exquisite size-dependent correlation factors in the 
Mott critical regime. 54 ' 
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Appendix: Variational Monte Carlo Method 

We briefly note the setup and condition of VMC cal- 
culations carried out in this paper. 

Because we need to optimize variational parameters 
up to the maximum number L, we use a correlated- 
measurement or optimization- VMC technique. 19 ' In the 
non-linear minimization process of energy expectation 
values, we adopt a quasi-Newton method, in which gra- 
dient vectors are effectively calculated, especially for 
bosons, by recently proposed formulae, 55 ' and Hessian 
matrices are approximately given by Broyden-Flecher- 
Goldfarb-Shanno formula, 56 ' the use of which does 
not affect the accuracy of optimization itself. In cod- 
ing, we refer to an algorithm offered by Ibaraki and 
Fukushima. 57 ' For wave functions with a few parame- 
ters, we use a simple linear optimization together. 

In both algorithms, parameters as well as energy con- 
verges typically after first several rounds of iteration 
with different fixed sample sets; in each set we gener- 
ate 2.4-2.5 x 10 5 particle configurations with Metropo- 
lis algorithm. After this convergence, we continue ex- 
cess rounds (15-90 times) of iteration in the optimization 
process with successively renewed configuration sets. We 
determine the optimized values by averaging the data 
obtained in the excess rounds; in averaging, we exclude 
scattered data beyond the range of twice the standard de- 
viation. Thus, the optimal value is substantially an aver- 
age of more than several million samples. The variational 
energy and significant parameters [g and /(2) etc.] are 
determined with sufficient accuracy in most cases, but 
accurate determination of insignificant parameters [/(r) 
with r > 7 etc.] is difficult, because E depends on them 
only very slightly, in other words, particle configurations 
determining them appear extremely rarely. Anyway, such 
parameters have little influence on E and other quanti- 
ties. Physical quantities are calculated with 2.4-2.5 x 10 5 
renewed configurations generated by the optimized pa- 
rameter sets. 

Since in Mott critical regimes, the global minimum 
becomes more competitive with other minima as L in- 
creases, accurate energy minimization sometimes be- 
comes not easy, especially for the triangular lattice with 
QWF and for *° p ^ and #j. This is the cause of scattered 
data points near U c /t in some figures. 
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